Nonlinear transport of Bose-Einstein condensates through mesoscopic waveguides 
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We study the coherent flow of interacting Bose-condensed atoms in mesoscopic waveguide geome- 
tries. Analytical and numerical methods, based on the mean-field description of the condensate, 
are developed to study both stationary as well as time-dependent propagation processes. We apply 
these methods to the propagation of a condensate through an atomic quantum dot in a waveguide, 
discuss the nonlinear transmission spectrum and show that resonant transport is generally sup- 
pressed due to an interaction-induced bistability phenomenon. Finally, we establish a link between 
the nonlinear features of the transmission spectrum and the self-consistent quasi-bound states of 
the quantum dot. 
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I. INTRODUCTION 



The development of microscopic trapping potentials 
for ultracold atoms has lead to a number of fascinat- 
ing experiments probing the behaviour of Bose-Einstein 
condensates on mesoscopic length scales. Examples in- 
clude the realization of a Josephson weak link between 
two condensates in a double well potential the mea- 
surement of interference and phase coherence between 
two spatially separate condensates 0, [l] , as well as the 
diffraction of a condensate from a magnetic lattice [3]. 
A convenient setup for such experiments is provided by 
"atom chips" Q where microscopic confinement poten- 
tials are created with the magnetic field that is induced 
by current-carrying electric wires mounted on top of the 
chip surface. This technique does not only allow one 
to produce microtraps, but also to create waveguide ge- 
ometries for cold atoms that can be rather flexible, and 
thereby opens the way to explore transport properties 
of cold atomic gases. Early experiments on atom chips 
did indeed focus on the propagation of a Bose-Einstein 
condensate along such a magnetic waveguide, where the 
condensate was transported in a controlled way by means 
of time-dependent magnetic fields Q or accelerated along 
the guide by means of a field gradient • 

The possibility to create such waveguides for cold 
atoms have stimulated a number of theoretical investi- 
gations on the transport physics of interacting matter 
waves, with particular emphasis on possible analogies 
with mesoscopic phenomena in the electronic context. 
This started with the attempt to define an atomic ana- 
log of Landauer's quantization of the conductance Q, 
and was continued by the generalization of the "Coulomb 
blockade" phenomenon to cold bosonic atoms propagat- 
ing through a quantum-dot-like potential M- More 
recent studies, which are based on an elaborate frame- 
work for the description of scattering processes of Bose- 
Einstein condensates (to be described in this article), in- 
clude the nonlinear resonant tr ansp ort of a condensate 
through atomic quantum dots [ill [l^ , the manifesta- 



tion or absence of Anderson localization in the trans- 
port through disorder potentials [l^, [T3l. as well as the 
transport of solitons through disorder [l^|. For their ex- 
perimental realization, these transport processes would 
require a coherent quasi-stationary flow of Bose-Einstein 
condensed atoms in the waveguide, which was recently 
realized in the context of optical guides 16] using the 
principle of "atom lasers" [rn |. 

From the theoretical point of view, the main compli- 
cation in the description of a quasi-stationary scattering 
process of a Bose-Einstein condensate obviously comes 
from the presence of the atom-atom interaction. In lead- 
ing order, the effect of this interaction is included in a 
nonlinear term in the Schrodinger-like Gross-Pitaevskii 
equation for the condensate wavefunction. In presence 
of a waveguide potential, providing a harmonic confine- 
ment in two (transverse) spatial dimensions and permit- 
ting free motion along the third (longitudinal) dimen- 
sion, an adiabatic treatment of the transverse degrees of 
freedom allows one to describe the evolution of the con- 
densate by means of an effective one-dimensional Gross- 
Pitaevskii equation as long as the confinement of the 
waveguide is sufficiently strong (such that the condition 
for the "ID mean-field regime" is satisfied [IBl)- This 
one-dimcnsional nonlinear wave equation does permit 
stationary solutions corresponding to condensates that 
propagate with finite velocity along the axis of the guide 
[ig |. As was shown by Leboeuf and Pavloff, these solu- 
tions can then be used in order to construct scattering 
wavefunctions of the condensate (with the appropriate 
outgoing boundary condition) in presence of finite-range 
perturbation potentials in the waveguide [20| . 

In contrast to the linear Schrodinger equation, the 
knowledge of stationary scattering states alone does not 
necessarily permit the prediction of the outcome of a 
given propagation experiment with Bose-Einstein con- 
densates. This is not only the case for the propaga- 
tion of finite wave packets (which obviosly cannot be de- 
composed into individual scattering eigenstates, due to 
the absence of the superposition principle in the Gross- 
Pitaevskii equation), but applies also to adiabatic in- 
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jection processes as performed in Ref. |l6|, where the 
waveguide is gradually filled with matter waves. Clearly, 
if such an adiabatic process leads to a quasi-stationary 
flow of the condensate (which actually need not be the 
case, as we pointed out in Ref. _13]), the corresponding 
scattering state necessarily satisfies the stationary Gross- 
Pitaevskii equation. However, not every scattering eigen- 
state of this nonlinear Schrodinger equation can eventu- 
ally be populated in this way: due to the nonlinearity, the 
eigenstates in the waveguide can be dynamically unstable, 
which means that they would disintegrate in the course of 
time evolution as a consequence of small deviations. Such 
dynamical stability properties cannot easily be inferred 
from the stationary Gross-Pitaevskii equation. Another 
nontrivial problem is, as we shall explain below, the de- 
termination of the incident flux of atoms that is associ- 
ated with a given stationary scattering state. This infor- 
mation is required in order to establish the connection 
to a given propagation experiment (where the incident 
current is typically under much better control than the 
net current during the propagation) and to determine the 
transmission coefficient of the scattering state. 

In view of these complications, it seems advisable 
to study waveguide scattering of Bose-Einstein conden- 
sates within the framework of the time- dependent Gross- 
Pitaevskii equation. While the straightforward numerical 
simulation of wave packet propagation processes is hardly 
feasible in the limit of spatially broad and energetically 
narrow wave packets (which would be required, e.g., 
for studying the energy-resolved transmission through 
atomic quantum dots), it is possible to directly simu- 
late the quasi-stationary injection process from an ex- 
ternal reservoir of Bose-Einstein condensed atoms into 
the waveguide, as it was experimentally performed in 
Ref. [TgI . Assuming that this reservoir is sufficiently large 
such that the effect of the back-action from the waveguide 
can be neglected, the dynamics in the waveguide is effec- 
tively described by an inhomogeneous Gross-Pitaevskii 
equation, which contains a source term that models the 
input of matter waves from the reservoir. This inho- 
mogeneous Schrodinger-like equation can be efficiently 
integrated with standard finite-difference methods, using 
absorbing boundary conditions in order to avoid artificial 
backreflections from the ends of the numerical grid. Typ- 
ically one would start with vanishing condensate density 
in the guide, and then time-integrate the equation while 
adiabatically increasing the source amplitude from zero 
up to a given maximal value. Glearly, this approach is 
rather close to the realistic experiment. By construction, 
it automatically yields, at the end of the propagation, 
scattering states that are dynamically stable (provided 
the flow remains quasi-stationary during the integration) , 
and it allows in a natural way to determine the transmis- 
sion of those states. We have successfully applied this 
approach to the transport of Bose-Einstein condensates 
through quantum-dot-like double barrier potentials [llj 
and through one-dimensional disorder potentials [l3| . 

The present paper is devoted to the detailed descrip- 



tion of this time-dependent approach to nonlinear waveg- 
uide scattering of a Bose-Einstein condensate, and to its 
relation with the existence of stationary scattering states 
of the condensate. To this end we briefly review in Sec.HIl 
the so called ID mean-field regime, set up the theoretical 
framework to study transport and scattering processes, 
and introduce concepts that allow to define transmis- 
sion and refiection coefficients for stationary scattering 
states that are solutions of a nonlinear wave equation. 
In Sec. Ill C[ the numerical method that is based on in- 
tegrating the time-dependent Gross-Pitaevskii equation 
in presence of a source term is explained. As a first ap- 
plication, the transmission spectrum of the condensate 
fiow through a quantum point contact consisting of a 
single potential barrier in the waveguide is discussed in 
Sec. Ill DI In Sec, mil we investigate the transport through 
a symmetric double barrier potential and we show in 
nil Al that the transmission spectrum exhibits an inter- 
action induced suppression of resonant transport. Fi- 
nally, in Sec. IIIIBI we develop an analytical description 
of the transport problem through the double barrier po- 
tential in terms of internal quasi-bound states. This es- 
tablishes a clear link between the nonlinear signatures of 
the transmission spectrum and the self-consistent quasi- 
bound states of the quantum dot. 



II. MEAN-FIELD APPROACH TO TRANSPORT 
OF CONDENSATES 

In the following we consider a coherent beam of Bose- 
Einstein condensed atoms at zero temperature, propagat- 
ing through a cylindrical waveguide with a finite-range 
scattering potential, given, e.g., by a constriction act- 
ing as a barrier potential for the beam. One of the 
aims of this work is to develop new methods to describe 
such propagation processes based on the Gross-Pitaevskii 
mean- field theory [21I, [22| . The mean-field dynamics of a 
dilute condensate can be described in terms of a macro- 
scopic order parameter, the condensate wave function 
^(r, t) , which obeys the nonlinear Gross-Pitaevskii equa- 
tion [231 



ih^^^ir,t) 



2m 



A + y(f) + t/o|*(f,i)|2 



*(f,t). 



(1) 

Low-energy scattering processes between two atoms in 
the condensate are predominately described by the con- 
tribution from s-wave scattering and lead to the non- 
linear term Uo\'^ {r, t)\'^ . Here Uq = ^ivfi^as/m is the 
interaction strength which is determined by the s-wave 
scattering length as and the mass m of the condensed 
bosons. The term V{r) in Eq. ([T|) is the external trap- 
ping potential experienced by the atoms. For the sake of 
definiteness we consider the experimentally relevant case 
of a condensate in a cylindrical harmonic waveguide with 
an additional scattering potential that is induced along 
the guide. Let x be the coordinate along the axis of the 
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guide and r = \J -V the cyhndrical radius associated 
with the transverse coordinates, then we assume V(f) to 
be of the form 



1 



V(r) = -TOw^r^ + V\\{x) 



(2) 



Here, the first term on the right-hand side is the trans- 
verse harmonic confinement of the guide with trapping 
frequency w and V||(a:) is the scattering potential paral- 
lel to the axis of the guide. could, e.g., consist of a 
single barrier that acts as a constriction for the conden- 
sate flow. Such a barrier can, for instance, be induced by 
irradiating a strongly focused blue-detuned laser beam 
onto the waveguide. 



A. ID mean-field regime 

In this subsection we derive an effective one- 
dimensional version of the Gross-Pitaevskii equation 
which is particularly suited to describe condensates 
in elongated waveguide structures. To this end, we 
adopt the adiabatic approximation method outlined in 
Refs. [H [l^, [13], where the condensate wave function 
can be cast into the form 



^'(a;, r) — 'ijj{x, t)(j){r, n). 



(3) 



Here, </> is the equilibrium ground state wave function for 
the transverse motion, normalized to unity 



1, 



(4) 



ip^x, t) describes the longitudinal motion, and the density 
per unit of longitudinal length is given by 



n{x, t) 



(5) 



We remark that this adiabatic ansatz involves a local den- 
sity approximation, in the sense that one assumes that 
the transverse motion depends solely on the local conden- 
sate density n{x,t) at position x. It was pointed out in 
Ref. [2^ that this approximation is justified if the trans- 
verse scale of the density variation is much smaller than 
the longitudinal one. This regime is certainly reached 
when the scale of variation of the longitudinal potential 
Vj|(x) is considerably larger than the harmonic oscillator 
length a± = ^Yij (um) of the radial transverse confine- 
ment. 

Inserting the ansatz ^ into the Gross-Pitaevskii equa- 
tion IT]) yields 
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r dr J 
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^muj'^r'^ + Uon{x,t)\(j)\'^ 



(6) 



We can identify the term in the square brackets as the ef- 
fective Hamiltonian Ht for the transverse degree of free- 
dom, acting on the wave function (f>, 



Ht4> — ^{n)4i. 



(7) 



The energy e(n) associated with the transverse state 
depends parametrically on the longitudinal density n. 
Thus, we obtain a pair of equations, one for the trans- 
verse, and one for the longitudinal dynamics of the con- 
densate. 



e(n)(j) 
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mu; r 



2m dx^ 



V\i{x)+e{n{x,t)) 



(8) 
(9) 



Eq. ([9]) is an effective one dimensional wave equation for 
the longitudinal order parameter which is particularly 
suited to describe a condensate in non-uniform waveg- 
uides. This regime is often denoted as the ID mean-field 
regime [Tsj. 

It remains to determine e(n). In the following, we as- 
sume that (p is the energetic ground state of Ht- In the 
so-called low density limit, a^n <C 1, the nonlinear term 
Uon\(l)\'^ in Eq. ([8]) is a small perturbation and a first- 
order perturbative solution of Eq. ^ yields 

e(n) = eo + C^o(</' ll^oPI 0o> ^ f-a + 2hhjasn, (10) 

where eo = fiw is the eigenenergy of the ground state 4>q of 
the unperturbed transverse Hamiltonian (eq is a constant 
energy shift which we drop in the following) . In the oppo- 
site large density limit, OsU ^ 1, the kinetic energy term 
in Eq. © can be neglected, and the so called Thomas- 
Fer mi ap proximation holds for the transverse wave func- 
tion [221, yielding 



1 



UoVn 



v/e(n)-14(r) 



(11) 



By imposing the normalization condition (|4]) to the 
Thomas-Fermi wave function (|lip we find in the high- 
density regime 



e{n) — 2hu}^ 



(12) 



At this point we remark that the validity of the Gross- 
Pitaevskii equation is restricted to the dilute gas regime, 
where the 3D density n^d fulfills nsaa^ < 1 [Ml, [ll]. 
This condition reads in the ID mean-field regime nag <IC 
{a_L/asf'/'^ (1/ = 1 in the low density regime and = 1/2 
for high densities [3). Typically a±/as is of the order 
10'^. This condition will be considered as always fulfilled, 
even in the regime of high longitudinal densities, when 
nug 3> 1. On the other hand, the weakly interacting ID 
Bose gas picture also breaks down at very low densities. 
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dimensional waveguide with a scattering potential in the 
ID mean-field regime. Starting point of our considera- 
tions is the effectively one-dimensional Gross-Pitaevskii 
equation To determine its steady solutions, we write 
ip{x,t) = A{x) ex.p[iS{x)]ex.p{~ifit), where A{x) and 
S{x) are real valued functions. The longitudinal density 
is n — A^, /i is the chemical potential of the condensate 
and V = {h/m){dS/dx) its local velocity. From Eq. © 
we obtain flux conservation n{x)v{x) = jt — const, and 
an equation of motion for the amplitude A(x) of the wave 
function 



m jf 



2m 



A" + - :^A + V» {x)A + e{n)A. (17) 



FIG. 1: Transverse energy e(n) (in units of hu) as a function 
of Usn (in units of fi? /m). The numerical result (solid line) 
coincides for large values a^n very well with the Thomas- 
Fermi result. The interpolation result agrees excellently with 
the numerical result for small values of OsTI and converges 
towards the Thomas-Fermi result for large asU. The inset 
zooms into the region of small a^n. The straight dashed line 
displays the perturbative result (|10|) . 



in the Tonks- Girardeau regime (see e.g. Refs. j25l. |26|. |27|. 
I2I]). This occurs in the regime nas <^ {as/o,±)'^ ~ 10~^ 
which we therefore discard from our present study. 

At the end of this section, we derive an analytical ex- 
pression that allows to interpolate e{n) between the two 
opposite limits na^ <^ 1 and na^ 3> 1. To this end we 
consider the ansatz 

e(n) =[a + /3{asn) + -fia^nf] ^'^ . (13) 

To determine the coefficients a, /? and 7, we expand 
Eq. (jl3p in the limit flgri ^ 1 to first order in a^n, 



e(n) = + ia-3/4^(a,n). 



{a,n) < 1. (14) 



In the limit a^n ^ 1 we keep only the quadratic term 
[asuf in Eq. E 



e(n) = 



(a^n) » 1. 



(15) 



The comparison of Eqs. (|14ll5p with Eqs. (|10ll2p yields 
a = S^w^, fi = Sfi.^'aj"' and 7 = 16?i^w'*, and the interpo- 
lation formula (fT3t reads 



j(n) = ^jK^uj"^ ^ A,W-uP-{asn). 



(16) 



This result can be compared with numerically computed 
values for e(n) [2^]. Indeed, as displayed in Fig.[Tl we find 
a good agreement between the interpolation result and 
the numerically computed values for the whole range of 
values of Ogn in between the two opposite density limits. 



B. Scattering states in waveguides 

In this subsection we study the stationary transport 
modes of a coherent condensate flow through a quasi one- 



In the following, we assume that the longitudinal poten- 
tial 1^1 (x) vanishes asymptotically in the "upstream" re- 
gion, i.e. for X — *■ —00, and in the "downstream" region, 
for X — > 4-00: y\\(x — > ±00) = 0. In accordance with this 
terminology, we consider an incident beam of condensate 
that propagates from x — > —00 to a: — > -t-oo (i.e., from 
the upstream to the downstream region) . 

In order to properly define the scattering problem, we 
first study the asymptotic behavior of the flow far away 
from the constriction, where V||(x) = 0. In this region, 
Eq. p?)) can be integrated once, yielding the first order 
equation of motion 



E = ^{AT 
Zm 



m Jt 
2A2 



+ liA^-£{n), (18) 



r 

with S{n) = / e{n)dh, 
Jo 



(19) 



where E is an integration constant. It was pointed out in 
Ref. [1^ that Eq. admits a simple interpretation in 
terms of classical dynamics, since it describes the energy 
conservation of a fictitious classical particle with "posi- 
tion" A and "time" x moving in the effective potential 



W{n) = (m jt)/{2n) + /in - S{n) 



(20) 



and the integration constant E corresponds to the total 
energy of the particle. Eq. is therefore integrable 
by quadrature (see Ref. [13] for a discussion of the low 
density regime o^n <C 1). 

The left panel of Fig. [5] displays the potential W{n) 
in the low-density regime where we have £{n) = 
with the effective interaction parameter g = 2fujjas'n'^. 
W{n) has qualitatively the same form in the high-density 
regime as well where £ is given by £{n) = 2gn^/2/3 
with g = 2hu}^fal . For weak and moderate coupling 
constants g (respectively g), W(ti) exhibits a local min- 
imum Emin — W{n\) at a low density n\ and a local 
maximum Emax = Win-i) at a high density ni. These 
extrema, at which the fictitious particle would be at rest 
forever, correspond to solutions of Eq. pS]) with con- 
stant density. They represent plane waves of the form 
ijj^ix, t) = y/n^exp{iki^x — i^t/h) (y = 1, 2) whose wave 
numbers are implicitly determined through the disper- 
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FIG. 2; Plot of the function W{n), here displayed for the low 
density regime £{n) = gr? 12, in the left panel. For given 
/i, j't and p, a beam of uniform density has either a density n\ 
(supersonic solution) or a density (subsonic solution). At 
a given classical energy E (with Emin < E < Emax, to assure 
bounded density oscillations ), n_ and are the minimum 
and maximum values of the cnoidal density oscillations, dis- 
played in the right panel. Energy values E close to (but lower 
than) Emax correspond to gray solitons. 



sion relation of the Gross-Pitaevskii equation , namely 

(21) 



2 nl 

as expressed in terms of the total current jt = hk^n,^/m. 
The solutions ipi and "02 are termed "supersonic" and 
"subsonic" , respectively, since the beam velocity is larger 
than the speed of sound of the condensate for tpi and 
smaller than the speed of sound for ?/'2 [^^l- The trans- 
port of particles at theses two solutions is dominated by 
the kinetic energy in the supersonic case, and by the in- 
teraction between the atoms in the subsonic case. We 
note that in the noninteracting limit, where e(n) is inde- 
pendent of n, the subsonic density n2 diverges and W(n) 
has only one finite extremum at the density ni. 



Solutions of Eq. ^ with E„ 



< E <Er, 



exhibit 



periodic density oscillations and correspond to a bounded 
motion of the fictitious classical particle. They are im- 
plicitly given through the integration of Eq. (fT51) . i.e. 



X — Xo 



A{x} 



y/2/m h dA 



, (22) 



where the amplitude A(xo) at the position xq determines 
the initial value for the solution of the differential equa- 
tion (fT7|) . For e{n) — gn/2 it was shown that the so- 
lutions of Eq. 122]) can be expressed in terms of Jacobi- 
EUiptic functions p^ . 

For our purpose, a qualitative characterization of the 
free solutions of the Gross-Pitaevskii equation is suffi- 
cient: Small deviations from the constant density value 
Til, e.g. small values of i? — W{ni) correspond to 
small sinusoidal density oscillations. Energy values close 
to (but lower than) the limiting classical energy value 
Emax — W(n2) correspond to gray solitons. In the in- 
termediate regime, between the limiting cases of small 
sinusoidal oscillations and gray solitons, the condensate 
density exhibits cnoidal oscillations. Energy values larger 



than -Emax lead to an infinite density at finite x and can- 
not be interpreted as physically meaningful steady-state 
solutions. We also note that the flat-density solutions co- 
incide, ni =712, when the potential W{n) exhibits a sad- 
dle point configuration. For the potential W{n) displayed 
in Fig. [21 such a saddle point configuration would, e.g., 
be encountered by increasing g while fi and jt are kept 
fixed. In the low density limit where e(n) = 2hujj_asn, the 
criterion for the existence of a saddle point configuration 
reads = 27'mj'^g'^; in the high density limit where 
e{n) — 2huj±y/asn, we find /i^ = 5^mj^g''/2^. Beyond 
these limits no stationary solutions exists any more. 

Finding stationary scattering states in presence of a 
finite scattering potential requires now to match two 
asymptotic density modes, each characterized by a sepa- 
rate integration constant E, in the upstream respectively 
downstream region. From general arguments on the dis- 
persion relation of elementary excitations of the Gross- 
Pitaevskii equation follows that the physically meaning- 
ful boundary condition for the steady-state solutions of 
Eq. I|17p demands a constant downstream density profile 
[l9|. The asymptotic downstream density should there- 
fore correspond either to ni or n2. In the present study, 
we intend to investigate the crossover from a noninter- 
acting to a weakly or moderately interacting system. We 
therefore focus on the regime of rather small conden- 
sate densities, respectively weak atom-atom interactions 
{nus ^ 1 and e(n) = gn); hence, the low-density down- 
stream solution rii will be relevant in the following. The 
high-density solution n2 exhibits qualitatively different 
features, such as solitonic transmission modes, and has 
been discussed in Ref. [20| . 

In analogy with the scattering problem in a non- 
interacting system we define a stationary scattering state 
as a solution of Eq. ^ of the form 



■0(x, t) = ip{x) exp{~i^JLt/h), 



(23) 



satisfying, in the downstream region, outgoing bound- 
ary conditions of the form ^{x) — y^ni exp(ifca:), with 
fc > given by ki as defined above. In order to deter- 
mine the scattering states for a given barrier potential 
Vj|(x), which vanishes at a; — + ±oo, and for given values 
for the total current flow jt and the chemical potential 
/i, we integrate the equation of motion p7p from the 
downstream to the upstream region with the "asymp- 
totic condition" A = ^Jn\ and A' = in the downstream 
region. This allows us to compute the density profile 
in the whole waveguide, and by computing the phase 
via S' [x) — mjtA'^{x)/h we determine unambiguously 
the stationary scattering state ipix). This procedure de- 
scribes the scattering process in terms of a so-called fixed 
output problem, because the outgoing current jt in the 
downstream region enters as a parameter in the asymp- 
totic boundary conditions that determine the scattering 
state Hill. 

There is only a small number of potential configura- 
tions, such as the square well or delta-peak barriers, for 
which the integration can be carried out analytically [12l | . 
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For the general case, it is convenient to rewrite Eq. p7l 
in terms of Hamilton-like equations of motion 



dp 
OH 
"dA 



P 



A^ 



2[A.-V||(x)-e(n)] A(24) 



where we introduced the canonical momentum p = 
(Y? jm^A! . These equations of motions can be deduced 
from the classical Hamiltonian 

n{A,p) ^^p^ + !^ + [^-V^^ ix)]A^ £{A'). (25) 

In the picture of the fictitious classical particle, l^|(a;) 
plays the role of a driving force which drives the par- 
ticle away from the minimum of the classical potential 
W{n). The classical energy, which is Ed = W{ni) in the 
downstream region, is altered by the amount 



AE 



+ 00 



Vii{x)A{x)A'{x)dx, 



(26) 



which yields the new classical energy value _E„ = Ed+AE 
for the upstream region. The asymptotic behavior of the 
scattering state is then fully determined by £'„ and Ed- 
The energy transfer AE is a measure for the amplitude 
of the density oscillations in the upstream region, i.e. 
increasing values of AE imply a larger backreflection. 

Our purpose is now to determine the reflection and 
transmission coefhcients, T and R, associated with the 
stationary scattering states. These quantities are natu- 
rally given by T = jt/ji and R=j^/ji, where ji, jt, and 
jr respectively denote the incident, transmitted, and re- 
flected current of the condensate. The determination of 
ji and jr, however, is a nontrivial task, since we can not 
simply decompose the upstream wave function into an 
incident and reflected plane wave component due to the 
nonlinearity of the Gross-Pitaevskii equation which does 
not permit the application of the superposition principle. 
We show now how the incident and reflected currents can 
nevertheless be deflned and calculated in a meaningful 
way. 

First, we briefly recall a method that has been sug- 
gested in Ref. [20| and was successfully applied in 
Ref. (isj . It allows one to determine approximate val- 
ues for T and R in the regime of small backreflections 
or small nonlinearities, by means of an approximate de- 
composition of the upstream density into an incident and 
reflected beam. We consider here the low density regime 
agU <C 1, e.g. e(n) = gn. In the upstream region, 
n{x) = A^{x) obeys the equation [see Eq. p^ ] 



with 



W{n) 



2n 



/in 



:9n 



(27) 



(28) 



We write the density in the form n{x) = ni+Sn(x), where 
dn{x) represents the density oscillations originating from 
back-reflections. Inserting this ansatz into Eq. ([27|l and 
introducing a new effective wave number 



K = fcv/1 - l/(2^2fc2) 



(29) 



(here, ^ = h/ ^2mnig is the condensate's healing length 
in the downstream region) and the characteristic scale 
5ni — m[Eu — W{ni)]/{ti^K'^) for density oscillations, we 
obtain 



dSn 
dx 



+ iK^Sn^ ^ SK^Smim + Sn) + ^^Sn^ (30) 



as an equation of motion for Sn{x). 

Until now, no approximation has been made. In the 
regime of small back-reflections, where |(5n|/ni <C 1 
holds, or small interaction parameters g (both limits are 
covered by the condition \Sn\/ni <^ k'^S,^, see Ref. '2^), 
we neglect the cubic term in Eq. pop . Thus, the equation 
of motion ([50)1 corresponds to the dynamics of a shifted 
harmonic oscillator, and its solution is given by 



n{x) = rti + Sni + \j2n\8n\ + (bn\)^ <zo^(2kx + 6*), (31) 



where Q is an arbitrary phaser. The density profile (j3ip 
is equivalent to that of the two countcrpropagating plane 
waves with wave vektor k 



ipi{x) 



Sni 

ni H — — exp(zra). 



Sni 



4>r{x) = \l CXp( — ZKX + i0), (32) 



yielding 



n{x) = \'4)i{x) + -tprix)]" 



(33) 



Identifying '0i(a;) as an incident and iprix) as a reflected 
wave component allows one to determine the transmis- 
sion and reflection coefficients through 

The approximate nature of Eq. becomes evident 
if we consider the conservation of currents. Computing 
the incident and reflected current components in the up- 
stream region yields 

/ Sni\ Hk . Sni fiK 

{ni + — ] — , Jr^^ , (35) 

\ 2/771 2 m 

whereas we find from the asymptotic downstream behav- 
ior of the wave function the transmitted current compo- 
nent jt — nihk/m. It is easy to see that the relation 
jt + jr — ji is exactly fulfilled only in the case of van- 
ishing atom-atom interactions, i.e. k = k. In the regime 
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FIG. 3: (color online) Adiabatic transition of the interac- 
tion parameter g for a proper definition of transmission co- 
efficients: The upper part of the figure displays the adia- 
batic variation of the position-dependent parameter g{x) from 
g = up to a maximal value g. The gray-shaded transition 
region between x\ and X2 in which g varies with position is 
assumed to be much larger than the typical periodicity of 
the condensate density oscillations. The lower part shows the 
density of a stationary scattering state in presence of a po- 
tential barrier, which is constant in the downstream region 
and displays oscillations in the upstream region (the position 
of the barrier potential is marked by the vertical line) . The 
nonlinear cnoidal oscillation of n between the barrier and X2 
is adiabatically conveyed, in the transition region between x\ 
and X2, into a sinusoidal oscillation in the interaction-free do- 
main on the left-hand side of x\. There, the wave function 
can be linearly decomposed into an incident and a reflected 
component. 



of weak interactions deviations from the current conser- 
vation are of the order ©[(fc^)^^] and the approximate 
approach becomes inappropriate for strong interactions 
or large backreflections. 

In order to overcome this problem, we consider a 
waveguide configuration in which the interaction strength 
g tends to zero for x — s- — oo and reaches a finite constant 
value in the region where the barrier potential is located 
(see Fig. [3]). We furthermore assume that the typical 
length scale on which g varies is much larger than the 
periodicity of the density oscillations. Such a variation 
of g can e.g. be achieved by decreasing the transverse 
confinement frequency lo of the waveguide or by tuning 
the scattering length as via a Feshbach resonance. 

Using once more the analogy with the dynamics of a 
classical particle, we introduce the effective "pseudo ac- 
tion" 



J = (p pdA = — 
m 



XQ+Ax 



[A'{x)]^dx 



(36) 



that is integrated over one spatial period Aa; of the 
upstream density oscillation (which would be given by 
Ax = Tr/k in the absence of the interaction). By use 
of Eq. the pseudo action can also be written in the 
form 

J^^\- ! ^ ^[Eu - W{n)]/n dn, (37) 
where , («+) is the minimal (maximal) density value of 



the oscillating upstream density. It is determined via the 
relation W{n±) — E^. Due to the theorem of adiabatic 
invariants, remains approximately constant along the 
waveguide as long as g is sufficiently slowly varied. It 
can, under this condition, therefore be evaluated at any 
position X, in particular also in the far-upstream region 
at X < xi where we have 5 = 0. There we can decompose 
the wave function in an incident and reflected part as 



ipix) = (ae*('=^+'^) + /3e-*'=^)e*'^ 



(38) 



with k — y/2mfl/h, where the amplitudes a, (3 and the 
phases ip, cf) are real. The wave function's amplitude reads 



A{x) = V"^ + + 2a(3 cos(2fca; -I- ^) , (39) 

and the canonical momentum p is given by 

, , , 2a[3ksin(2kx + w) 

p{x) = —A'{x) = ^ ^' . (40) 

TO + /32 + 2a/3 cos(2fcx ^ ip) 

By use of dA — A'dx, we evaluate Eq. ((36|l as 

J <j) pdA,^2f3^h'^kT:/m. (41) 

Using the fact that the incident and reflected currents 
read ji = hka^/m and jr = hkp'^/m in the far-upstream 
region, we can obtain the reflection and transmission co- 
efficients via 



R 



T 



Jr 
ji 



(3^ 



l-R=l 



J 



1 



J 



2TThjt 



(42) 



Eq. (|42p unambiguously assigns a reflection and a trans- 
mission value to each scattering state that is a solution 
of the nonlinear wave equation ^ . This deflnition repre- 
sents a natural extension of the concept of transmission 
for nonlinear scattering problems. 

For the practical computation of the transmission value 
associated with a given scattering state, it is sufficient 
to evaluate the integral (I37|) numerically in the near- 
upstream region (i.e. for x ~ a;2 in Fig. [3]) where the 
extremal densities n± can be found by solving W{n±) = 
Eu- This means that the adiabatic variation of g does 
not need to be included at all in the calculation; it is suf- 
ficient to take into account a short spatial domain in the 
upstream region within which the condensate exhibits 
a couple of density oscillations. Computing the trans- 
mission by use of Eq. circumvents the approximate 
character of the relation ([M)) and is therefore also valid 
in the regime of strong atom-atom interactions as well 
as for large back-reflections. In Fig. [4] we compare the 
approximate with the "exact" expression for the trans- 
mission, determined by Eqs. ([M)) and ([1^ respectively, 
for a condensate with a moderate nonlinearity that en- 
counters a potential barrier in the guide. For small large 
back-reflections both results coincide, whereas for large 
back-reflections the approximate formula (|34p systemat- 
ically overestimates the transmission. 
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FIG. 4: (color online) Comparison between the approximate 
and exact transmission values, calculated with Eqs. (I42|l and 
(|34|) respectively, for a moderately interacting condensate that 
propagates through a potential barrier, with the dimension- 
less parameters = 3, g = 1/2, jt — 1, and with vari- 
able barrier height. The transmission is plotted as a func- 
tion of the classical energy transfer A_E, which is a mea- 
sure for the back-reflection in the upstream region. For in- 
creasing back-reflections, the approximate result (|34|) overes- 
timates the "true" value of the transmission coefficient given 
by Eq. (|42|l . 



C. Time-dependent transport processes 

So far, we restricted our considerations to station- 
ary scattering solutions of the time-independent Gross- 
Pitaevskii equation. A severe problem is the fact that the 
mere existence of a stationary scattering state does not 
imply that this state is dynamically stable and can be 
populated in a time-dependent scattering process. This 
is not only true for the propagation of a finite wave packet 
(which obviously cannot be evolved by an expansion 
in terms of stationary solutions of the Gross-Pitaevskii 
equation , due to the absence of superposition principle), 
but also concerns the limiting case of a quasi-stationary 
flow that is generated by an adiabatic injection of the 
condensate into the waveguide. This affects, as we shall 
discuss later on, the resonant transport of a condensate 
through a double barrier potential, where the dynami- 
cal stability properties of the scattering states become 
crucial for their population. 

In view of this complication we now describe a method 
based on the time- dependent Gross-Pitaevskii equation , 
which allows us to simulate a realistic propagation pro- 
cess. This equation is integrated in presence of an inho- 
mogeneous source term, located at a position x = xq in 
the upstream region and emitting monochromatic mat- 
ter waves. The source term simulates the coupling of 
the waveguide to a large reservoir of a Bose-condensed 
matter at a given chemical potential /x, from which mat- 
ter waves are injected into the waveguide (see Fig. [5]). 
The effective nonlinear wave equation that governs the 
time evolution of the condensate wave function ip{x, t) is 
therefore given by 



ih 



d4>{x, t) 
dt 



I^^^V {x) 
2m dx"^ ' 



mx,t)[^ 



-\-S{t) 5{x — Xq) exp (— i/ii/ft). 



7/'(x,t), 



(43) 




scattering potential V(x) 



FIG. 5: (color online) A reservoir of Bose-condensed matter 
with a given chemical potential ^ is locally coupled at the 
position X(i to a waveguide with a scattering potential. The 
reservoir emits a plane matter wave in both directions into the 
guide. Hence, a coherent beam, with current ji, propagates 
towards the barriers of the potential where the condensate is 
partially reflected, with the current jV, and partially trans- 
mitted, with the current jt. 



where the time-dependent coupling strength between 
the waveguide and the reservoir is contained within the 
source amplitude S{t). The interaction parameter g need 
not be constant, but may be considered to position- 
dependent as well, in order, e.g., to simulate the adia- 
batic transition from a noninteracting to an interacting 
guide as depicted in Fig. [3l In this work, we restrict our- 
selves to the case where g is constant in the vicinity of 
the finite-range scattering potential. 

Before studying time-dependent scattering processes 
in a waveguide with a finite scattering potential, it 
is instructive to consider first stationary solutions of 
Eq. (j43| for the particular case of a homogeneous waveg- 
uide, i.e. V||(a::) = 0, and a constant source amplitude 
S{t) = So- In this case, there exist plane wave solutions 
tp{x,t) = ip{x)e~^^'^ with constant density n — \ip{x,t)\'^. 
To demonstrate this, we switch to the Fourier space 
by introducing the Fourier transformed wave function 
ij{q,t) = J exp{iqx)^{x,t)dq. Then, Eq. takes the 
form 



ih 



d_ 

"dt 



2m 



gn ijj{q,t) = 5*0 6" 



tqxo 



-ifj,t/h 



This equation admits solutions of the form 
i>{q,t) 



9m 9^ 



(44) 



(45) 



Here, we introduced the wave vector k via the relation 
h?k'^ — 2m(^—gn). By transforming back to the position 
space, we find solutions where the source term emits in 
both directions the monochromatic wave 



(46) 



with the wave number k being self-consistently defined 
by 



2ml ^- g 



(47) 
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FIG. 6: Illustration of the relation (|50p between the source 
amplitude So and the densities ni,2 of a homogeneous flow. 
The lower branch corresponds to the supersonic solutions with 
density ni and the upper branch to the subsonic solutions 
with density n2 . The insets illustrate the configuration of the 
classical potential W{n) for the different types of solutions. 



The density that is associated with the wave function 
reads 



\So\^m 



(48) 



Evaluating the quantum mechanical current operator 
shows that the source emits the current 



Ji 



± 



\So\ 



(49) 



with " + " for X > xq and "— " for x < xq. Inserting 
Eq. into Eq. immediately yields the plane-wave 
dispersion relation /i = mjf / (271^) + gn. 

As already discussed in Sec. IIIBI (see Eq. ((2T|) ) this 
equation admits two flat-density solutions, n — ni and 
n = 712, corresponding to a supersonic and a subsonic 



propagation of the condensate, respectively. 
Eq. (gSl) in the form 

1 



"1,2 = ^ (mt vV 



2gm\So\yh' 



Rewriting 



(50) 



allows one to compute the two densities ni_2 that are pos- 
sible for a given value of the source amplitude 5*0. This 
relation is illustrated in Fig. [6l the lower branch con- 
tains the supersonic solutions and the upper branch the 
subsonic solutions. The value Smax = hfi/y^2gm corre- 
sponds to the saddle point configuration of the classical 
potential W{n); for source amplitudes larger than this 
threshold, no stationary solutions of Eq. (^5]) are pos- 
sible. In the limit of noninteracting particles, g = 0, 
only the supersonic branch survives (because the speed 
of sound is zero) and Eq. (|50p takes the simple form 
n = \So\^m/{2h'^^i). 



Now we study the time evolution of iIj{x, t) in pres- 
ence of a variable source amplitude S{t). Here, the sce- 
nario of an initially empty waveguide that is gradually 
filled with matter waves is of peculiar interest as this cor- 
responds to the experimentally realistic situation where 
the condensate is initially confined in a microtrap (play- 
ing the role of the reservoir) and then smoothly released 
to propagate into the waveguide. To simulate such a 
process, we propagate ipi^i t) by numerically integrating 
the wave equation (j43p in presence of an adiabatic in- 
crease of the source amplitude S{t) from S{t = 0) = up 
to a given maximal value Sq, with the initial condition 
ip{x, t = 0) = 0. The amplitude S{t) is increased adia- 
batically in order to ensure that, at any instant during 
the propagation, the wave function in the guide remains 
as close as possible to a stationary scattering state of the 
form ijj{x) exp^—ifit/h). Quantitatively this means that 
the typical time scale AT on which the amplitude S{t) 
increases is much larger than the characteristic time scale 
T = h/fi that is associated with the chemical potential 
fi of the source: AT ^ r. As we are studying an in- 
finitely extended scattering problem, we have to impose 
absorbing boundary conditions in order to avoid artificial 
back-reflection at the boundaries of the numerical grid. 
Details on these absorbing boundaries, which are taken 
from Ref. fs^ and adapted to account also for a finite 
nonlinearity, as well as on the numerical integration pro- 
cedure are given in Appendix |B1 

In a first step we discuss the filling of the waveguide 
in absence of a scattering potential, i.e. for Vj|(a;) = 0. 
Fig. Elja-c) displays the time-evolution of the wave func- 
tion ip{x,t) by a series of snapshots showing the density 
at different times. For the sake of definiteness we chose 
S{t) = 5o[l — cxp(— i/AT)] , which provides a smooth 
evolution towards the desired final value S{t — > oo) = 5*0. 
We find that for propagation times t ^ AT the calcula- 
tion converges towards the flat density (Fig.[71;) that cor- 
responds to the stationary plane wave (|46)) at the source 
amplitude S — Sq. The bottom part of Fig. [7] shows the 
real and imaginary parts of the wavefunction -0 at the 
time t — AT during the filling process. These panels 
clearly illustrates that the source emits a plane-wave like 
solution of the form A{x,t) exjp[—i^t + ik{x,t)x] where 
A{x,t) and k{x,t) vary slowly with position and time. 

It is instructive to display the evolution of the conden- 
sate density as a function of the time-dependent source 
amplitude S. Fig. [8] shows this evolution for different 
values of AT (dashed lines). We notice that these curves 
approach the lower branch of the relation ([50]) (solid line 
in Fig. [51 see also Fig. [5]) if we reach the limit AT ^ r. 
We therefore deduce that the adiabatic filling of an ini- 
tially condensate-free waveguide can only populate sta- 
tionary solutions that correspond to a supersonic flow; 
hence, the flnal condensate density is given by n = ni 
as deflned by Eq. ([50)) . In analogy to the fixed output 
problem discussed in Sec. Ill Bl the implementation of the 
source term therefore allows one to investigate the trans- 
port of the condensate in terms of a so-called fixed input 
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(a) 



(b) 



(c) 






FIG. 7: Time evolution of the condensate density during the 
adiabatic increase of the source amplitude (the source is lo- 
cated at the vertical dashed lines). The panels (a-c) show 
three snapshots of the condensate in a waveguide without 
scattering potential: at (a) t = 0.1 AT, (b) t = AT, and (c) 
t — 10 AT. The bottom part of the figure shows the real and 
imaginary part of the wave function whose density is displayed 
in panel (b). The panels (d-f) illustrate the scattering of the 
matter waves at a repulsive barrier potential (gray-shaded re- 
gion). Panel (f) clearly shows that a stationary scattering 
state is populated in the long-time limit t ^ AT. 



problem, where the incident current ji that is emitted 
into the guide parametrizes the process. The fixed in- 
put approach is much closer to experimental situations 
because the current that is injected into the guide is typ- 
ically under much better control than the total transmit- 
ted current jt- 

In a second step, we consider a scattering process 
in presence of a barrier potential V||(a;). Due to the 
partial backscattering of the condensate at the bar- 
rier, the dynamics becomes more complex as compared 
to the potential-free case. Nevertheless, for weak or 
moderate nonlinearities the wave function ip(x,t) is 
found to converge towards a stationary scattering state 
ipix) exp(— during the adiabatic increase of the 
source amplitude towards its final value 5*0, as illustrated 
in Fig. [Tljd - f). During the gradual filling of the guide, 
the condensate is partially reflected at the barrier, which 
leads to the oscillating density pattern between the bar- 
rier and the position of the source in the upstream re- 
gion. On the right-hand side of the barrier, in the down- 
stream region, the density is flat in the long-time limit 
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FIG. 8: (color online) Evolution of the condensate density as 
a function of the source amplitude S for three different values 
of the time scale AT in which S is ramped to its maximal 
value 5*0 (dashed and dashed-dotted curves). For an increas- 
ing ratio AT/t with r = h/fi, the curves converge towards 
the supersonic branch of Eq. ((50} (solid line). For S — > So 
the supersonic scattering state with constant density n = ni 
is reached. 



t ^ AT, which reflects the fact that the wave function 
Tpix) is given there by an outgoing plane wave of the 
form = ^exp(iA:a;). We checked that the state V'(a^) 
that is reached at the end of the propagation fulfills the 
stationary Gross-Pitaevskii equation, i.e., the wave func- 
tion's amphtude A{x) = \ip{x)\ is a solution of Eq. pT|) . 

Once we populate a stationary state, we have another 
straightforward access to the transmission coefficient T 
in the nonlinear scattering problem: T is given by the ra- 
tio of the transmitted current jt, evaluated through the 
current operator in the downstream region, to the current 
that would propagate through the waveguide in absence 
of the barrier potential, which is the current ji that is di- 
rectly emitted from the source. This approach provides 
another natural extension of the definition of transmis- 
sion coefficients to nonlinear wave equations. Hence, the 
numerical method introduced in this section allows not 
only to calculate scattering states that are dynamically 
stable and can be populated in a realistic propagation 
process [s^, but provides also a straightforward access 
to transmission coefficients for a fixed input problem. 

We point out that in the nonlinear case, convergence 
towards a stationary scattering state is not always guar- 
anteed. Indeed, studying the transport of condensates 
through a waveguide with an extended disorder region by 
means of the method described in this section revealed 
that, beyond a critical interaction strength respectively 
a critical length of the disorder region, the transport 
process generally remains time-dependent and stationary 
states are not populated [l3| . 



D. Transport through a quantum point contact 

As a first and simple example we study the transport 
of a Bose-Einstein condensate through a quantum point 
contact. We consider a waveguide with a constriction 
given by a single repulsive Gaussian barrier potential 
V||(a::) = Vq exp(— x^/cr^) which can be experimentally 
implemented by focusing a blue detuned laser beam in 



11 



its transverse ground mode onto the waveguide [3J| . For 
the sake of definiteness we consider in the following a con- 
densate of ^''Rb atoms (to = 1.45 x lO^^^kg, Ug = 5.77 
nm) flowing through a waveguide with transverse trap- 
ping frequency w = 27r x 10'^ s~^ that corresponds to 
a harmonic oscillator length aj_ = 0.34 fim. It is con- 
venient to measure energies in units of fuv, lengths in 
units of aj_ and particle currents in units of w. In these 
units the interaction parameter reads g = 0.0,iAhuja±. 
For the longitudinal extension of the barrier we assume 
a — 2a^ ~ 0.7/iTO (which would be at the limit of experi- 
mental realizability) and its height is chosen as Vq = 3 fiw. 

In a first step, we investigate the transport process in 
terms of a fixed output problem: We calculate scatter- 
ing states by integrating the stationary Gross-Pitaevskii 
equation, 



2m dx^ 



2 I 2 

-X j a 



^g\i>{xt 



V^(x),(51) 



for given values of the chemical potential /i and the trans- 
mitted current jt from the downstream to the upstream 
region (where a supersonic density n(x — > oo) — rii is as- 
sumed in the downstream region). Eq. allows us to 
compute the corresponding transmission coefficient from 
which we can deduce the incident current via ji = jt/T. 
Varying the transmitted current allows us to compute 
the jt — ji current characteristics which is displayed for 
/i = 3fiw in the inset of Fig. [9] For noninteracting parti- 
cles {g = 0) the jt — ji characteristics is linear, because 
the transmission coefficient T does not depend on the 
particle current, whereas for non- vanishing interaction 
parameters the jt — ji characteristics shows a nonlinear 
behavior and displays increasing deviations from the lin- 
ear case with increasing particle currents. This means 
that the presence of repulsive interactions suppresses the 
transmission through the quantum point contact with 
growing current. 

It is now easy to switch from the fixed output to the 
fixed input problem where the incident current ji is kept 
constant. To this end, we basically have to invert the 
jt — ji characteristics, in order to determine the total 
current jt and the corresponding transmission coefficient 
T that result from a given incident current ji. This can 
be done in a unique way in the present case, since the 
jt — ji characteristics is monotonous and therefore allows 
one to unambiguously assign to each value of jt a unique 
incident current ji. Computing the current characteris- 
tics for different values of fji allows one then to obtain the 
transmission spectrum, i.e. the transmission coefficient 
T as a function of the chemical potential /i at a fixed 
incident current ji. Transmission spectra of the point 
contact for different incident currents are displayed in 
Fig. m Qualitatively, we find that in presence of repulsive 
interactions the spectra resemble strongly the spectrum 
for a single particle: for chemical potentials considerably 
smaller than Vq the transmission tends to zero, whereas 
for fjL much larger than Vq we reach a regime of perfect 
transmission. In the intermediate regime, we clearly see 
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FIG. 9: (color online) Transmission spectrum of the conden- 
sate flow through a quantum point contact for different values 
of the incident current ji (/i in units of fioj). The solid lines are 
found by evaluating the stationary Gross-Pitaevskii equation, 
the values denoted by crosses are obtained by integrating the 
time-dependent Gross-Pitaevskii equation in presence of the 
source term (the dashed line displays the result for a noninter- 
acting condensate) . The inset shows the jt — ji current char- 
acteristics for a condensate flow without interactions (g — 0, 
dashed line) and in presence of interactions {g — 0.034:huj, 
solid line), at /i = Shu. 



that increasing particle currents ji yield a moderate sup- 
pression of the condensate flow through the point con- 
tact. This is attributed to the fact that the presence of 
the repulsive interaction leads, at fixed /i, to a reduction 
of the available kinetic energy, which in turn reduces the 
probability for the atoms to penetrate the barrier. 



So far, the computation of the transmission spectra 
was based on the stationary Gross-Pitaevskii equation. 
As a complementary access, we apply the method based 
on integrating the time-dependent Gross-Pitaevskii equa- 
tion with source term. For each value of fj,, the wave func- 
tion was propagated according to Eq. ([ii]) in the pres- 
ence of an adiabatic increase of the source amplitude S 
up to the maximum value Sq that corresponds to a given 
incident current ji. For the considered range of incident 
currents ji we find stationary scattering states at the end 
of the propagation. As shown in Fig. [9l the results for 
the transmission obtained from the time-dependent in- 
tegration (marked by blue crosses) coincide with the re- 
sult based on evaluating the stationary Gross-Pitaevskii 
equation. Hence we can conclude that a gradual filling 
of the guide populates precisely those scattering states 
that are eigenmodes of the stationary problem, and that 
these stationary states are dynamically stable. 
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III. TRANSPORT THROUGH A DOUBLE 
BARRIER POTENTIAL 

Now we study the particularly interesting propagation 
process of a Bose-Einstein condensate through a sym- 
metric repulsive double barrier potential which can be 
seen as a Fabry-Perot interferometer for matter waves. 
This setup was first discussed by Carusotto and La Rocca 
[l3| who proposed to use a combination of optical 
lattices for the realization of this bosonic quantum dot. 
In the context of atom chips, a double barrier potential 
could also be implemented by suitable geometries of mi- 
crofabricated wires on a multilayer chip geometry [ssj . 
Another straightforward implementation relies on two 
blue-detuned parallel laser beams, crossing transversely 
the waveguide. Assuming the laser beams to be in the 
lowest transverse mode, this setup creates a potential ge- 
ometry with two Gaussian shaped barriers. For the sake 
of definiteness, we consider this latter case and assume a 
double barrier given by 



(52) 



Here, a is the width of one barrier and L is the distance 
between the barriers. 

For a flow of noninteracting particles it is well known 
that the transmission spectrum of a symmetric double 
barrier potential exhibits Breit- Wigner resonances (36| 
which are related to resonant transport states. In our 
context, these resonant states can be defined as station- 
ary scattering states of the condensate (see Eq. (^5]) ) 
that exhibit perfect transmission. In the following, we 
investigate to which extent resonant transport through 
such a double barrier potential can be achieved for an 
interacting condensate, and how interactions modify the 
transmission spectrum. 



A. Resonant transmission spectra 

We now compute transmission spectra for the double 
barrier potential (|52p by applying the same methods that 
have been employed to find the spectra of the quantum 
point contact in Sec. Ill D[ using again the same units 
that were already introduced there. In the following, we 
consider a condensate with effective interaction strength 
g (which will be varied to investigate the effect of an in- 
creasing nonlinearity) , a waveguide with transverse trap- 
ping frequency w = 27r x 10^s~^, and a double barrier 
potential (j52[) with the parameters Vq — 1.1 fno^ cr = a^, 
and L = 4.25a_L. We study the transport of the con- 
densate in terms of a fixed input problem, with incident 
current ji = 1.0 cj. The influence of the atom-atom in- 
teraction on the transmission spectrum is exemplarily in- 
vestigated in the vicinity of the energetically lowest res- 
onance which has one density maximum in between the 
two barriers (see inset in Fig. fTU]) . In Ref . |ll| we showed 



that qualitatively similar results are also found for higher 
resonances. 

First we compute transmission spectra by use of 
the integration method based on the stationary Gross- 
Piatevskii equation, respectively Eq. (jl7p : The spectrum 
is, as in Sec. IIIDl determined by calculating stationary 
scattering states for given jt and /i, and the incident cur- 
rent of the scattering states is computed via Eq. (|42p . 
Finding the value of jt that results from a given ji is an 
optimization problem that can be solved systematically 
by analyzing the jt — ji current characteristics. In the 
linear case, 5 = 0, we obtain a Breit- Wigner resonance 
a.t fj, — 0.389huj corresponding to the energetically lowest 
resonance state fFig. [TU|. 

Now we consider the case of a weak atom-atom in- 
teraction, g — 0.002?iu;a^. As the most striking result, 
we find, close to the resonance, a multivalued transmis- 
sion spectrum where two further solutions appear for 
0.419 < n/ihv) < 0.472. These solutions join together 
to form a resonance peak that is asymmetrically dis- 
torted towards higher values of the chemical potential 
[sj. The resonant state, which is found at /i = O.A72huj, 
coexists with a low-transmission state, as depicted in 
the central panel of Fig. [TOl The asymmetric distortion 
becomes even more pronounced for an increasing inter- 
action strength. This is shown in the bottom panel of 
Fig. [TO] where we dipslay the spectrum in the vicinity of 
the first resonance for g = Q.01huj±. 

It is instructive to trace the evolution of the jt — ji 
characteristics in the vicinity of the onset of the mul- 
tivalued subzone in the spectrum. In contrast to the 
monotonously increasing current characteristics that we 
found for the quantum point contact (Fig. [9]), the char- 
acteristics of the double barrier potential shows a more 
complex behavior, where it is not always possible to un- 
ambiguously attribute to each incident current ji one sin- 
gle transmitted current jt. Fig. [TT] shows that for values 
of fj, below the critical chemical potential from which on 
three branches coexist, the current characteristics inter- 
sects only once the horizontal line that represents the 
fixed incident current ji = 1.0 w. Above this critical value 
of fi three intersection points are found, corresponding to 
the three coexisting scattering states. 

Our findings are characteristic for a bistability phe- 
nomenon, similar to processes in nonlinear optics (ssj 
and in the electronic transport through quantum wells 
(39I l40j . It is crucial to know which branches of the 
transmission spectrum are actually populated in a real- 
istic experimental situation in order to decide if resonant 
transport is possible in presence of a finite interaction 
strength. To this end, we recalculate the transmission 
spectrum with the time- dependent integration approach, 
which simulates, at given value of /i, the adiabatic release 
of the condensate from the reservoir into the waveguide. 
As explained in Sec. IIIDl this method provides another 
straightforward access to the transmission values, and 
stationary states that are selected by this method auto- 
matically satisfy the criterion that they are dynamically 
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FIG. 11: (color online) Current characteristics for the double 
barrier potential in the vicinity of the onset of the multival- 
ued subzone of the spectrum (for g = 0.002 }uja±). Below 
the critical fi — OAlQhuj, the ji — jt characteristics (dashed- 
dotted line) intersects only once the horizontal line indicat- 
ing the fixed incident current ji = 1.0 a;. Above this critical 
value, three intersection points are found (dashed line). At 
fi — OAlQhij, the current characteristics (solid line) exhibits a 
tangent to the horizontal line. The intersection points marked 
with filled (green) points correspond to scattering states that 
are populated during a time-dependent propagation process. 



FIG. 10: (color online) Transmission spectra of the double 
barrier potential ai g — (upper panel), g = 0.002 hija± 
(middle panel), and g = 0.01 hu!a± (lower panel). The (green) 
solid lines show the transmissions of all scattering states, cal- 
culated by the "stationary" method based on Eq. (|17p . that 
exist at the incident current ji = 1.0 w of the matter- wave 
beam. The dashed lines display the spectra obtained from 
the time-dependent integration approach. The inset shows 
the longitudinal atom densities (in units of a~[^) of the first 
resonant state and the coexisting low-transmission state for 
g = 0.002 huja±_ (the position x is given in units of a_i_). The 
gray-shaded curves indicate the positions of the two barriers. 
The (red) dots (marked by the arrows) designate the posi- 
tions of the resonant state and the low-transmission state in 
the transmission spectrum. 



stable and can be populated in a realistic propagation 
process. 

The dashed lines in Fig. [10] show the result of this cal- 
culation. While a perfect agreement with the method 
based on the stationary Gross-Pitaevskii equation is 
found for 3 = 0, the time-dependent approach repro- 
duces, for g 7^ 0, only the lowest branches of the spectra 
in the multivalued region. This apparently implies that 
the asymmetrically distorted peak structure is essentially 
inaccessible in the propagation process that is considered 
here. We therefore conclude that resonant transport, 
which would necessarily require the population of such 
a distorted peak, will generally be suppressed in pres- 
ence of finite interactions, and only the low branches of 
the spectrum which have rather low transmission will be 
populated. Qualitatively, this behavior of the nonlinear 
system can be understood by comparing the "internal" 
interaction energy (evaluated within the internal region 



of the double barrier) 



-L/2 



Eint = g I dx 

-L/2 



(53) 



of the resonant with the one of the coexisting low- 
transmission state. The system can minimize -Ejnt by 
realizing a state with a low particle density in between 
the barriers. As displayed in the inset of Fig. (TUl this 
favors the low-transmission state. 

To conclude this section, we remark that a tempo- 
rary enhancement of the transmission of matter waves 
near the resonance can be achieved by a variation of 
the external potential during the propagation process. 
In Ref. [ll| we devised a temporal modulation scheme 
where the potential is shifted with time according to 
V{x) V{x,t) = V{x) - Vb(t). Specifically, such a 
modulation can be induced by illuminating the scatter- 
ing region with a red-detuned laser pulse, where Vo{t) > 
would be determined by the detuning and the intensity of 
the laser. In the case of an adiabatic modulation of V , the 
wave function 'ip{x, t) remains, at each time t, close to the 
instantaneous scattering state that is associated with the 
external potential V{x, t) — or, equivalently formulated, 
close to the scattering state for the potential V{x) at the 
shifted chemical potential /i + Vb {t) . As soon as + Vb {t) 
is raised above the critical chemical potential from which 
on the transmission spectrum becomes multivalued, the 
wave function follows continuously the upper branch of 
the resonance and evolves into a near-resonant scatter- 
ing state with high transmission. This state turns out 
to be dynamical unstable, and the wave function decays 
after a typical lifetime of the order of several milliseconds 
towards a low-transmission state 11 1. 
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B. Transmission in terms of quasi-bound states 

In this subsection, we present analytical and numer- 
ical evidence that the distortion of the resonance peak 
arises indeed due to the nonlinearity-induced level shift 
of the self-consistent quasi-bound state within the atomic 
quantum dot. We describe, for this purpose, our system 
in a similar way as in the well-known scattering matrix 
approach [ilj, namely by a discrete "bound" (or quasi- 
bound) state within the quantum dot that is weakly cou- 
pled to two symmetric continua of unbound "lead" states 
in the up- and downstream regions of the waveguide. In 
contrast to the situations for which the scattering ma- 
trix formalism was originally developed [4l| , we consider 
here nonlinear dynamics within the quantum dot, which 
is described by the Gross-Pitaevskii equation. As was 
pointed out above the outcome of a given scattering pro- 
cess is, in this case, not completely independent of the 
"history" of the process, i.e., of the way in which the con- 
densate is injected into the waveguide. Different scatter- 
ing states might, specifically, be populated if the chemical 
potential is adiabatically varied in different ways during 
the propagation [ll|. To account for this complication, 
we formulate our nonlinear scattering theory in a time- 
dependent way, namely by considering the asymptotic 
propagation of a spatially broad (and energetically nar- 
row) wave packet that is injected onto the quantum dot 
from the left (upstream) lead. The population of the 
wave packet that exits the scattering region in the right 
lead gives naturally rise to the transmission coefficient. 

As starting point, we subdivde of the Hilbert space 
Ti. into a subspace Tip containing discrete bound states 
within the quantum dot region, and two other subspaces 
Hl/u containing continuous states in the left and right 
leads of the waveguide. This subdivision can be formally 
achieved by means of the Feshbach projection method 
p3 |. where those subspaces are defined by the projec- 
tion operators Pl = 9{xl — x), Pr — 9{x — xr), and 
Q = 1 — Pl — Pr- Here xl and xr are suitably chosen 
positions that mark the left and right boundaries of the 
quantum dot, and 9 denotes the Heavyside step function. 
As an essential ingredient of the Feshbach formalism, dif- 
ferent boundary conditions (i.e., of Dirichlet or Neumann 
type) are imposed within and outside the dot, which al- 
lows one then to shift the boundary contributions from 
matrix elements of the Laplace operator to appropriate 
sides of the spatial cuts at x — xl/r, in such a way that 
the operator T of the kinetic energy remains Hermitean 
within each subspace, but exhibits finite coupling matrix 
elements across the boundaries (see, e.g., Ref. [i^ for 
more details). Choosing Dirichlet boundary conditions 
within the resonator and Neumann boundary conditions 
in the leads, these matrix elements would read 

ii'Rim = ^rR{xR)<p'ixR,) (54) 

(Vl|T|0) = -^^1(xl)4>'(xl) (55) 
2m 



for wave functions (/>(a;), ipLix), and iPr{x) defined within 
the subspaces Ho, Hl, and Hr, respectively. Without 
loss of generality, we set x^ = —a and xr = a in the 
following, where a = L/2 denotes the position of the 
maximal barrier height. 

We now make the assumption that the nonlinearity 
can be neglected in the lead regions outside the quantum 
dot, which should be valid at weak interaction strengths 
and which is motivated by the fact that close to reso- 
nance the density within the double barrier potential is 
strongly enhanced as compared to the leads. We further- 
more assume that only one quasi-bound state, namely 
the local "ground state" of the quantum dot, apprecia- 
bly contributes to the scattering process, which is in- 
deed the case in our specific double barrier potential ([52]) 
where "excited" quasi-bound states are energetically lo- 
cated above the barrier height. Neglecting the contribu- 
tion of those excited states, we make the ansatz 

/>oo 

ij{x,t) = / dEA'^{t)(j)'l,{x) + B{t)Mx) 

JQ 

/>oo 

+ / dEA^{t)cf>^{x) (56) 
Jo 

for the wave function, where (f>Q G Hq denotes the above 
quasi-bound state and (f)^ G ^l/_r are the energy- 
normalized continuum eigenstates within the left and 
right lead, respectively, at energy E. Inserting this ansatz 
into the Gross-Pitaevskii equation yields the equations 

z;i^4/^(i) = EA'^J''{t) + VEB{t) (57) 

,n±Bit) = ^,o{\Bit)\^)Bit) 

+ dEVE[Al:{t)+A§{t)] {58) 
Jo 

for the amplitudes Aj^, A^, and B. Here, 

^o(|S(i)P)^Mr+.9|SWP (59) 

with 

~g = gf \M^)\^dx (60) 

J — a 

represents the population-dependent chemical potential 
of the quasi-bound state, and 

Ve ^ ^^'o{a)^^{a) = '^M-a)^Ei-a) (61) 

denotes the coupling matrix element between 0o and 
4>^^^- We assume here, without loss of generality, that 
the wave functions (j)o{x), (j)j^{x), and 4>e{x) are real- 
valued and that the continuum eigenfunctions exhibit the 
symmetry-related property 4>e{^) = 0b (~^)- 

As appropriate initial state for the quasi-stationary 
scattering process, we consider a spatially broad Gaus- 
sian wave packet that is injected from the left-hand side 
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onto the double barrier potential. This wave packet is 
explicitly written as 



ip{x, t^) = a exp 



2^ 



(62) 



with Xe = xq/c^ and Ce = cro/e^ for xo,aa > 0. Choosing 
the initial time in the asymptotic past according to 
te — ~mxe/{hk), the wave packet will, in the limit e — > 
0+, evolve into the plane wave 



(63) 



at finite times t, with the incident chemical potential /i = 
h'^k'^/{2m). Using the fact that the energy-normalized 
continuum eigenfunctions are, in the asymptotic spatial 
region a; 3> a, given by 



2m 

TTh?kf 



cos{kEX + ipE) (64) 



with kE = V2mE/h and with a potential-dependent 
phase we obtain the initial amplitudes 



Wkl 



-a exp 



-yAkE-kf 









X exp 


+iXe (^kE - 


+ i^PE 



(65) 



and B{t,) = A§{t,) = for e 0+. 

Equation (j57p can now be formally integrated yielding 

-^Ve f B{t')e-'''^*-''y^dt' . (66) 

Inserting this expression into Eq. ([58)) leads to the equa- 
tion 

ihj^Bit) ^ ,,o{\Bit)\^)Bit) 

/' dt' Blt')e-'^^'-'"^'''K{t - t') 

h .h 



POO 

+ / dEVEA%{Qe- 



/"OO 

K{t) = / dEV^e-'^'^-''''^/'' . 
Jo 



-iE(t-U)/h (-g-.^ 

for the bound component, with the Kernel 

(68) 



In the limit e — > 0, the last term on the right-hand side of 
Eq. (|F7|) is evaluated as 5'e~*^*/'' with the effective source 
amplitude 



S 



277 h?k 



(69) 



This suggests that the time-dependence of the bound am- 
plitude is, in the quasi-stationary case, dominated by the 
exponential factor e"*^*/''. 

This latter information permits now to evaluate the 
second term on the right-hand side of Eq. ((67)) : if 
B(t') exp{ifit' /h) varies much more slowly with time than 
K{t — t'), we can justify the approximation 



dt' Bit')e~"''-*-'"^/''K{t-t') ~ B{t) / drKir) 
t, Jo 



(70) 



where the energy shift (5^ and the rate 7^ are, respec- 
tively, given by the principal value integral 



and by the expression 



(71) 



(72) 



Omitting the small shift in the following, we obtain 
the equation 



rn^^Bit) 



f^oi\Bit)\')~-hj,\Bit) 



(73) 



for the bound component B{t), which exhibits strong 
analogies to a nonlinear damped oscillator model that 
is subject to a periodic driving. Obviously, stationary 
solutions of Eq. (|75)) are of the form 



B{t) = Boe-'"'/^ 



(74) 



where the bound amplitude Bq satisfies the self- 
consistent equation 



Bn = 



(75) 



For the noninteracting case g — 0, one can show that this 
solution is necessarily realized after a transient propaga- 
tion time of the order of 7^^- 

Inserting this stationary solution into the equation 
for the transmitted component finally yields 



-2m 



Y2^-ifj.{t-tc)/h 

-A^At.) (76) 



while A^{t) would, for E ^ n, vanish in the limit e ^ 0. 
We therefore obtain the transmission coefficient through 



r(M) 



\A^it)\' ^ 

\Aj:iuw [fi-f,oiW)f + {hj,/2y 



(77) 



In the noninteracting limit 17 — > 0, this expression de- 
scribes the Breit-Wigner profile of a single resonance 
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peak at ^ = Ho- Indeed, if the decay rate is sufBciently 
small around this resonance, we can safely approximate 
7m by 7mo in the relevant energy range \^ - fio\ < hj^.^. 
Then T(fi) is given by a perfect Lorentzian centered 
around fj, — fio with the width h'y^g. At finite g 0, 
however, T may exhibit several branches for a given value 
of /i, due to the implicit relation (|75p between the bound 
component Bq and the incident chemical potential fi. 

We now aim at reproducing the numerically calculated 
transmission spectrum (see Fig. [TU]) through Eqs. ((77|) 
and ((75|) using information that is obtained from the cor- 
responding decay problem [H, IH, 113, HI], namely 
the chemical potential and the instantaneous decay rate 
of the local quasi-bound state at given population |i?oP- 
The latter quantity can also be derived from Eq. ([67|) . 
now in absence of the incident wave Aj^{tg) and with the 
initial population -B(io) = ^o- Taking into account the 
fact that the dominant time-dependence of B(t) is, in 
this case, given by exp[—ifi{\Bo\'^)t/h] for not too long 
evolution times t, we obtain 



^n^^Bit) 



f^o{\B{t)\')--njo{\B{t)\^))B{t) (78) 



as equation for the bound component B{t), with 
7o = 7po(|-BP)- Clearly, Eq. (|78|l describes a non- 

exponential decay of the condensate in the quantum dot, 
which is explicitly given by the equation 



-^o[Nb{t)]Nt.{t) 



(79) 



where the decay rate varies adiabatically with the re- 
maining population Nb(t) = \B{t)\'^ of the quasi-bound 
state. Such nonexponential decay processes of Bose- 
Einstein condensates were discussed in detail in Refs. [i^, 
113, HI], where the instantaneous decay rates 70 (.^Vf,) at 
various populations Nj, were used to predict the time 
evolution of the quasi-bound population through the nu- 
merical integration of Eq. ([751) ■ 

In analogy with the noninteracting case, we now re- 
place 7p 7o in Eq. (|77|) . which approximately 
interpolates between the decay rate of the weakly popu- 
lated quasi-bound state at /i = Mo(0) (which is naturally 
given by 70 (0)) and the decay rate near maximum of the 
shifted resonance peak. Using this approximation, the 
equation for the transmission coefficient reads 



2 ' 



(80) 



where the quasi-bound population Nb implicitly depends, 
via Eqs. (|75)) and ([55)1 . on the incident chemical potential 
/i and the incident current ji = hk\a\'^/m according to 



Nbit) 



hMNb)/2 



(81) 



- ^^oiNb)V + [hjoiNb)/2Y 

As in the corresponding decay problem [H, [13, HI], 
we now need to know the instantaneous chemical po- 
tentials fio{Nb) and decay rates -joiNb) at given quasi- 
bound populations Nb in order to calculate solutions of 




FIG. 12: Chemical potential fio and decay rate 70 of the 
quasi-bound state within the double barrier potential, cal- 
culated as a function of Ntg with Nb the population of the 
quasi-bound state and g the effective one-dimensional inter- 
action strength. In practice, /lo and 70 were computed at 30 
equidistant values of Ntg within ^ Ntg ^ 1.5, and cubic 
interpolation was employed to obtain intermediate values of 
fio and 70 for the self-consistent solution of Eq. (fST)! . fj,o, hjo, 
and g/a are given in "natural" energy units of huj. 



this set of equations. We apply for this purpose a real- 
time propagation method which is based on the numer- 
ical integration of the "homogeneous" time-dependent 
Gross-Pitaevskii equation (i.e., without the inhomoge- 
neous source term) in presence of absorbing boundaries. 
Starting from an appropriate initial condensate wave 
function (which should approximate quite well the reso- 
nance state to be calculated) , and renormalizing the wave 
function after each propagation step to satisfy the con- 
dition 



\'iP{x)\'^dx = Nb 



(82) 



within the quantum dot, one indeed obtains, after a suf- 
ficiently long propagation time, convergence towards the 
lowest decaying state of the system. The scaling factor 
that is needed to perform the renormalization (|82p gives 
then rise to the decay rate 70 = 70 (Ah) of the quasi- 
bound state, while the chemical potential ^0 = fJ-oiNb) 
of the decaying state can be extracted from the expecta- 
tion value of the nonlinear Gross-Pitaevskii Hamiltonian. 
In practice, it is sufficient to compute /io and 70 in this 
way for the equidistant values Nbg = 0, 0.05, 0.1 .. . of the 
population A^^, and to use cubic interpolation in order to 
determine intermediate values of /iq and 70. 

With this information, the possible self-consistent val- 
ues of the quasi-bound population can be computed by 
applying a numerical root-search method to Eq. (j8ip at 
given chemical potential /i and given incident current ji. 
The resulting occupation numbers Nb are then inserted 
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FIG. 13: (color online) Transmission spectra of the double 
barrier potential ai g — (upper panel), g = 0.002 huia± 
(middle panel), and g — 0.01 hu;a± (lower panel). The solid 
line shows the transmissions of all scattering states, calcu- 
lated by the "stationary" method based on Eq. ((17} , that ex- 
ist at the incident current ji — lu of the matter-wave beam. 
The dashed line is obtained from self-consistent solutions of 
Eq. (|81|l at — lu, which are inserted in the expression (|80p 
for the nonlinear transmission coefficient. The good agree- 
ment confirms the one-to-one correspondence between quasi- 
bound states of the atomic quantum dot and resonance peaks 
in the transmission spectrum (/x in units of hu)) 



in the expression ([50)1 for the transmission coefRcient. As 
shown in Fig. [T21 a distorted resonance peak is then ob- 
tained for g > 0. Apart from a slight overestimation 
of the peak width, this peak agrees quite well with the 
peak structure that would be formed through the trans- 
mission coefficients of all possible stationary scattering 
states at the above incident density. This ultimately 
confirms the one-to-one correspondence between quasi- 
bound states of the atomic quantum dot and resonance 
peaks in the transmission spectrum. 

It is worthwhile to note that self-consistent solutions of 
the quasi-bound populations can also be found in a dif- 
ferent way, namely by iteratively inserting approximate 
expressions for Ni, into the right-hand side of Eq. ([5T|) 
starting with Ni, — 0. This approach would effec- 
tively mimic the quasi-stationary propagation of a Bose- 
Einstein condensate through the initially empty quantum 
dot. In agreement with the time-dependent propagation 
approach based on the inhomogeneous Gross-Pitaveskii 
equation (see Sec. Ill C[) . only the lowest branch of the 
distorted resonance peak is populated in this way. This 
again underlines that the framework used in this section 



is intrinsically suited to take into account time-dependent 
effects and might therefore be used to predict the out- 
come of specific propagation processes. 



IV. CONCLUSION 

We have presented analytical and numerical results for 
steady and time-dependent flows of repulsively interact- 
ing Bose condensed atoms through mesoscopic waveg- 
uide structures. To this end, we described a theoretical 
framework that is suitable to study transport and scat- 
tering processes in the ID mean-field regime. In this con- 
text we introduced a non-perturbative method to extend 
the concept of transmission and refiection coefficients to 
nonlinear wave equations. On the other hand, to pre- 
dict the behavior of the condensate fiow under realistic 
experimental conditions, it is necessary to study time- 
dependent transport processes. We developed for this 
purpose a numerical method based on integrating the 
time-dependent Gross-Pitaevskii equation in presence of 
a source term that simulates the coupling of the waveg- 
uide to a reservoir from which a quasi-stationary flow of 
condensate is smoothly released into the guide. 

The approach was flrst applied to the transport 
through a single quantum point contact, where we found 
as a main result that an increasing nonlinearity leads to a 
distinct reduction of the transmission. Much more com- 
plex behavior was found for the condensate flow through 
a double barrier potential. Here, the atom-atom interac- 
tion induces a bistability phenomenon of the transmitted 
flux in the vicinity of resonances, which manifests as a 
strong distortion of the transmission peaks. By means of 
the time-dependent integration scheme, we demonstrated 
that resonant transport will consequently be suppressed 
in a realistic propagation process. However, as we showed 
in Ref. (llj . a suitable variation of the external poten- 
tial during the propagation process can enhance the fiow 
to reach a near-resonant state on finite time scales. Fi- 
nally, an analytical description of the transport problem 
through the double barrier was developed, which estab- 
lishes a clear link between the nonlinear signatures of 
the transmission spectra and the properties of the self- 
consistent quasi-bound states of the quantum dot. Simi- 
lar results were recently obtained in Ref. [i^ as well. 

Our numerical approach based on the inhomoge- 
neous time-dependent Gross-Pitaevskii equation can be 
straightforwardly generalized to describe scttering pro- 
cesses in multidimensional geometries. It can certainly 
be applied also to more complex scattering potentials, 
involving more than two barriers. In that case, however, 
we do not expect that the calculation always converges 
towards a stationary scattering state, even if the source 
amplitude in the inhomogeneous Gross-Pitaevskii equa- 
tion is varied on a very long time scale. This was demon- 
strated in our study on the transport of Bose-Einstein 
condensates through one-dimensional disorder, where we 
found that randomly generated disorder potentials of fi- 
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nite range will generally give rise to permanently time- 
dependent scattering processes at finite interaction, as 
long as the length of the disorder region exceeds a critical 
interaction-dependent value [l^ [ij] . Interestingly, this 
cross-over between quasi-stationary and time-dependent 
scattering, arising for disorder samples with lengths be- 
low and above this critical value respectively, correlates 
with a transition from an exponential (Anderson-like) to 
an algebraic decrease of the average transmission with 
the sample length [3], which indicates that the deple- 
tion of the condensate during the propagation process 
might play a prominent role there. We note in this con- 
text that the effect of depletion can to a certain extent 
be accounted for within the framework of our approach, 
namely through the implementation of the microscopic 
quantum dynamics approach introduced by Kohler and 
Burnett (50| in combination with an external source [5l| . 

The results that were obtained in this work are related 
to other fields of nonlinear physics as well, such as nonlin- 
and the electronic transport through quan- 
[40| , where similar observations on resonant 
transport were made. In the context of Bose-Einstein 
condensates, the realization of a quasi-stationary flux 
of interacting matter waves though scattering potentials 
that are defined on microscopic length scales still rep- 
resents a formidable experimental challenge. There are, 
however, promising advances in this direction, such as 
the atom-laser-like injection of a condensate into an op- 
tical waveguide [l^ as well as the scattering of a sta- 
tionary condensate in presence of a moving obstacle [s^l ■ 
Such advances should, in combination with detection 
techniques for single atoms (ssl . [s^ (which would allow 
one to measure very low transmissions), make it possible 
to experimentally investigate the role of interaction in 
mesoscopic transport processes from a new perspective, 
namely the one of cold bosonic atoms. 
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(in the following we set for simplicity h = I, m = 1 
d 

with the effective nonlinear Hamiltonian 

H{x, f) = -i ^ + V{x) + g\^{x, t)\\ 



(Al) 



(A2) 



which we want to integrate for a given initial state 
'0(x,io) of the condensate. In order to compute the 
time evolution of the condensate wave function ip{x,t) 
for t > to, we subdivide the time interval t — Iq into n 
discrete time steps of the size At — {t — to)/n, and use 
an implicit Crank-Nicholson integration scheme [ssl ] to 
propagate the wavefunction from one time step to the 
next one. The effective time evolution operator U for 
one discrete time step At is then given by [s^ 



U{t + At,t) = 



1 



1 



jH{x,t)At 



1 



-H{x,t)At 



(A3) 



The representation (jA3[) of lA is unitary and thus con- 
serves the norm of the wave function tp. The implicit 
integration scheme for the wave function reads then 



1 + —H ] i}j{x,t+At) 



l-^i/)^(x,t). (A4) 



We expand the wave function on a discrete lattice with 
N lattice sites by introducing the grid basis 



X3 



1 : Xi 



\Ax ^x < Xj 



: otherwise. 



(A5) 



with Ax = {Xjyiin Xrnax)/^ ■ Here, Xmin and Xmax S-IG 

the boundaries of the finite grid. The wave function then 
reads 
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ip{x,t„) 



N 



(A6) 



where ijj" = i^{xj, t„) is value of the wave function at the 
position Xj of the j'th lattice site (the index n labels the 
discrete times, t„ = tQ + nAt). Using the finite-difference 
representation for the kinetic part of H{x,t), we find 



l±^i/)v(^„t„)^V^;±^x 



2Ax^ 



(A7) 



APPENDIX A 

In this appendix we describe the numerical integration 
procedure of the time-dependent Gross-Pitaevskii equa- 
tion and the implementation of the source term into this 
integration scheme. We consider the equation of motion 



with Vj = V{xj). By introducing t})'^ = (?/;". ..-(/'"...V'tv) ' 
the lattice representation of Eq. (jA4[) finally reads 



where we define 



Di = 



<^ = Da^^DiV^" , (AS) 



(A9) 





, D2 = 
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and the N x N matrix representation of D1.2 reads 



with 



±a 1 =F l3j ±a 



iAt iAt / I , 



2 V Aa;2 



+ y,+g|^;n.(AlO) 



4Ax2' 



Hence, the integration of Eq. (|Aip reduces to the solution 
of a system of hnear equations with a tridiagonal matrix. 

So far, our integration scheme uses the value of tp" 
at the beginning of the integration step. This neglects 
the fact that the effective Hamiltonian (|A2|1 is implicitly 
time-dependent due to the presence of the nonlinear term 
g\il^{x,t)\'^ . Thus, it would be appropriate to use a more 
precise estimate for this nonlinear term, which is some- 
how averaged over the timestep At leading from t„ to 
tn+i- This problem can be handled by using a predictor- 
corrector-like scheme which was already successfully ap- 
plied in [srj . In this scheme, each integration step is 
done twice: First, we propagate the wave function from 
time tn to time tn+i using ip^ in the nonlinear term, in 
order to obtain a predicted wave function ■0"+^. Then, 
we repeat this integration step but using now the aver- 
aged value ^[V'" + in the nonlinear term, yielding 
a corrected wave function ■0" . 

Now we consider the presence of the source term. The 
equation of motion reads therefore 

d 

'i-Q^^{x, t) ^ H{x, t)il;{x, t)+S{t) exp{-ifit) S{x). (All) 

Working with a grid representation of the wave function, 
it is convenient to approximate the (5-function by 

R{x) = -J- \e(x + Ax/2) - eix ~ Ax/2)] , (A12) 
Ax 

where is the Heavyside step function. Before including 
the source term to the finite difference scheme, we esti- 
mate the error that is introduced by this approximation. 
To this end, we study the steady-state solutions of the 
wave equation 



.dtp 



1 

2dx^ 



Tp = So R{x) , (A13) 



that are obtained in the limit t — s- cx). The Green func- 
tion that is associated with the stationary equivalent of 
Eq. (|A13|) is given by 



G{x - x') 



ik 



(A14) 




FIG. 14: (color online) Real and imaginary parts (black 
solid lines) of the steady-state plane-wave solution obtained 
by integrating the time-dependent Gross-Pitaevskii equation 
with the numerical source term (|A12|I using the time step 
At — h/ (50/x) and the grid spacing Ax = A/20 with A — 2n/k 
the wavelength of the condensate. An excellent agreement 
with the exact analytical result (|A13|I (red dashed lines) is 
found. The source is located at the position x = xq. 



with k = ^/2{p- gn) (see Sec. HTdI Eq. dH])). Hence, 
the ansatz 



+ 00 



dx' ^e*^-|"-"'l R{x') 
ik 



(A15) 



yields a solution ipj^{x) of Eq. (jAlSp . Evaluating this 
integral yields 



V'„(x) 



2^0 



-ikx 



s\n{kAx/2) : x < 



ik^Ax 



2 

1 - ^os(fc2:) : 1^1 < ^ 

At 

e*'=^sin(A:Aa;/2) ■ x > — 



(A16) 

which converges towards Eq. (|A14p in the limit Aa; 0. 
The result (jA16[) can serve as an estimate for the relative 
error that is done by approximating S{x) with R(x): 
we obtain 

2sin(fcAx/2) PAx^ kAx , 
^ = 1 r-^ — =i — r— if —:— < 1. (A17) 



kAa 
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The relative error therefore scales quadratically with the 
grid spacing Ax and becomes negligible for reasonably 
small values of Ax. 

The above considerations justify the implementation 
of the source term at the position xy through the dis- 
cretized form 



s^i = s{t„)eM-i^^tn)SJ,, 



(A18) 



where 5jj' = 1 if j = j' and otherwise. In the presence 
of the source term, Eq. (|A8p is modified and reads 



D2^"+^ + 5" = DiV^" ^ rP''+^ = D^^Dii/)'" - 6") 

where the components of the vector 6" are given by 
iAt 



(A19) 
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In Fig. [Ml we compare the exact result (|A13p to the nu- 
merically computed plane-wave solution that is obtained 
in the limit t ^ cx) by simulating the gradual filling of 
a waveguide without scattering potential, = 0. In- 

deed, we find an excellent agreement between the numer- 
ical result and the exact plane-wave solution (IA13|) if we 
choose, e.g., Ax = A/20 (with the wavelength A ~ 27r/fc) 
and At = ft/(50/i). 

It is worthwhile to mention that in the presence of 
strong nonlinearities (for values of g considerably larger 
than in this paper) and strong backreflection, a nonlin- 
ear back-action between the reflected matter wave and 
the source term can occur. As a consequence, the trans- 
mitted current depends not only on the source amplitude 
Sq but also on the position of the source. In such a situa- 
tion, it is advisable to implement the adiabatic transition 
scheme that is displayed in Fig. [3] where g vanishes in 
the far-upstream region. By positioning the source term 
there and by choosing a sufhciently large transition re- 
gion, one can avoid this nonlinear back-action and ensure 
that the wave function is adiabatically conveyed from a 
linear wave to a nonlinear scattering state obeying the 
Gross-Pitaevskii equation . 



APPENDIX B 

In the numerical treatment of time-dependent scatter- 
ing processes in open quantum systems, one often en- 
counters the problem of defining physically meaningful 
boundaries at the edges of the computational domain. 
The naive, straightforward expansion of the wave func- 
tion on a finite spatial grid generally leads to an artifi- 
cial backscattering of the wave function from the bound- 
aries of the grid, which makes it impossible to simulate 
infinitely extended scattering states. This problem can 
be circumvented by introducing complex absorbing po- 
tentials in the vicinity of the grid boundaries (see, e.g., 
Ref. [3) J which should be designed such that they ab- 
sorb the outgoing flux as best as possible without af- 
fecting the dynamics inside the scattering region. An 
alternative method, which was introduced by Shibata for 
the linear Schrodinger equation [s^ . consists in the def- 
inition of absorbing boundary conditions (ABC) at the 
edges of the grid, which are formulated in order to per- 
fectly match outgoing plane waves with a specified dis- 
persion relation. This method is particularly suited for 
quasi-stationary propagation processes where the out- 
going part of the wave function is well described by a 
plane monochromatic wave. We show in this Appendix 
how this approach can be numerically implemented, and 
how the effect of a moderate nonlinearity in the Gross- 
Pitaevskii equation can be taken into account. 

We first discuss the absorbing boundary conditions for 
the Schrodinger equation (with h — I and to = 1) 



1^ 

2&?2 





(2a^ ) 


1 









co-y 



FIG. 15: (color online) The positive branch of the dispersion 
relation of a plane wave (black line) is approximated by a 
linear function (straight blue line). The parameters a\,a2 
are chosen such that the wave numbers of the plane waves to 
be absorbed lie within the momentum interval Afc. 



where 14 is a constant potential which is independent 
of the position x. This equation admits plane-wave so- 
lutions ip{x,t) — ^e^*^'^*"'^^^ satisfying the dispersion 
relation 



(B2) 



The "-|~" and "— " branches of Eq. (|B2p correspond to 
plane waves that propagate to the right- and left-hand 
side, respectively. Thus, the ABC should satisfy the dis- 
persion relation given by the branch of Eq. (jB2[) 
at the right boundary and the " branch at the left 
boundary of the grid. 

We derive now so called "one-way wave equations" on 
the basis of the dispersion relation l|B2p . which we will 
implement at the boundaries of the grid and which locally 
allow for wave propagation only in the outgoing direction. 
To this end, we make use of the duality relations 



d_ 
dt 



d_ 

dx 



ik 



(B3) 



which is going to be inserted into the dispersion rela- 
tion (jB2[) . Unfortunately Eq. (jB2[) is nonlinear in /j, and 
cannot be straightforwardly converted into a linear differ- 
ential equation. To circumvent this problem, we approx- 
imate Eq. (jB2|) in the vicinity of the chemical potential 
of the wave to be absorbed by the linear function 



k = ±- 



Q!2 — ai 



/2ai a2\/2ai - aiV2a2 

M ± : r (B4) 



Oi2 — Oil 



(see Fig. [T5|) . The parameters Q;i,a2 are chosen such 
that Eq. (jB4p is a good approximation to the dispersion 
relation (|B2p within the interval Afc = \J2ai — \/2ol\ 
around the central wave number i(^/2a2 + \/2a7). By 
use of the duality relations (jB3p . Eq. (|B4[) is transformed 
into the one-way wave equation 



-oTI-F+^e W'(a:,t), (Bl) 



. '9'/' / . 1 9 gi 
at \ gi ox gi 



V', 



(B5) 



21 



with 



51 

92 



±- 



± 



a2 - ai 



a.2 — oil 



(B6) 



Implementing these one-way wave equations at the 
boundaries of the grid (see below) leads to a very good 
absorption of plane waves with wave numbers k satisfy- 
ing y/2ai < fc < \/2a2- In Ref. [s^l it was demonstrated 
that also wave packets of the form tp = Ai exp(ikjx) 
can be absorbed if all wave numbers in this superposition 
lie within the above interval. 

It is straightforward to see that the one-way equa- 
tions (|B5|1 absorb plane waves also in the presence of the 
nonlinear term gjV'P- This is evident for the special case 
of a constant density: ip(x,t) = y^exp(— i/ii/fi, ± ifcx) 
with the dispersion relation 



±y/2{^ - gn). 



(B7) 



is obviously a solution of the Gross-Pitaevskii equation. 
A comparison of Eq. (|B7p with Eq. (|B2p reveals that the 
term gn can be identified as a constant effective potential. 
Hence we set Ve = gn for a proper absorption of the plane 
wave. 

We now generalize this result for plane waves whose 
parameters are slowly varying in time and position. This 
case is of high relevance for our work since the gradual 
filling of the guide with matter waves leads to the popu- 
lation of a scattering state whose outgoing parts, which 
have to be absorbed at the boundaries of the grid, exhibit 
slowly varying amplitudes and phases. We consider 



(B8) 



where A(x, t) and S{x, t) represent the local amplitude 
and phase, respectively, of the wave function. Locally, at 
position X — xq^ we can expand the phase according to 

S{x, t) = 5(xo, t) + k{xo, t){x - xo) +0[{x~ a;o)^](B9) 

with k{xo, t) = dxS{x, t)\x=XQ In the limiting case where 
Aix^ t) and S{x, t) vary on time and length scales that are 
considerably larger than l//i and l/fc(a;o,t), respectively 
(for X ~ xq and for all times t), Eq. (jBSp locally takes the 
form of a plane wave with a slowly varying amplitude and 
wave number. Under this condition, we find at a given 
position Xo at any time the local dispersion relation 



k{xQ,t) = ±v/2(^ - gn{xo,t)). 



(BIO) 



withn(xo,t) = |A(a;o,t)p. Hence, fc(a;o, t) parametrically 
depends on t via the condensate density at the position 
xq which is supposed to be at the boundary of the grid. 
By adjusting the values of ai and a2 such that \/2ai < 
k{xo, t) < \J2a2 is satisfied for all times i, the wave V' is 
absorbed at the edge of the lattice. 



N-2 N-1 N 

• — e — • 

^ Ax ^ \ 
Xr-2Ax Xr-Ax 



Additional point x 



FIG. 16: Sketch of the right lattice boundary. The additional 
point at position x allows for a proper implementation of the 
absorbing boundary conditions in a grid representation of the 
wave function. 



We now outline how to incorporate the ABC into the 
lattice representation (|A6|) of the wave function. Here we 
consider exemplarily the right-hand side boundary Xr — 
xn of the grid, where the wave function has to obey 
Eq. (|B5|) with the upper {+) sign in the definition (|B6p 
of the prefactors. The idea is now to replace the equation 
for the boundary component i/'at of the state vector -0", 
i.e. the last component in the equation (|A8|) . by the finite- 
difference version of the one-way wave equation (jBSp . To 
this end, we need a finite-difference expression for the 
derivative dx4'ix,t)\x=xr at the grid boundary. Since Xr 
is the last grid point, this expression can only be obtained 
in an asymmetric way with respect to Xr, namely through 
the difference between ^{xr,t) and 7/1(2;^ — Ax,t). This 
would lead to the equation 

■^[i^i^r, t+At)- yj{Xr, t)] = (Ve - ^{Xr,t) 
i ^p{Xr,t) — ip{Xr — Ax,t) 
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Ax 



(Bll) 



which was also used in Ref. [32|. 

The asymmetric structure of Eq. (jBlip introduces a 
small but systematic error in the propagation of the wave 
function, since the value and the derivate of ^ are, strictly 
speaking, computed at different positions, namely at Xr 
and at the intermediate point x = Xr — Ax/2, respec- 
tively. This problem can be circumvented by replacing 
Eq. (jBlip with the analogous equation for the wave func- 
tion ip{x,t) evaluated at this intermediate point x (see 
Fig. [T6|) . There we have 



^^(x,i) 



l\)(Xr-, t) — -ipiXr — Ax, i) 

Ax 



(B12) 



as "exact" (i.e., symmetric) finite-difference expression 
for the derivative, and the value of the wave function at 
this additional point is obtained through 



1 



'lp{x, t) ~ - [^{Xr,t) + IpiXr — Ax, t)] 



(B13) 



Inserting these expressions (jB12p and (jB13l) into Eq. (|B5p 
leads to a symmetric finite-difference equation for 
ip{xr,t) and •ip{xr — Ax,t) where the value of the wave 
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function at the auxiliary point x does not explicitly ap- 
pear any longer. In the grid representation, this finite- 
difference equation reads 



2Ai 



in-"' 



rN - rN-i) 



giAx 
51/ 



m - rN-i) 

(B14) 



Eq. (jB14[l allows for a straightforward incorporation into 
the matrix representation (|A9[) : the modified matrices 
Di 2 read at the right-hand side edge of the numerical 
grid 



Di 



a 1 — (3n-2 a 

a 1 — Pn~i ol 

V 73 74/ 



( . \ 

-a 1 + I3n-2 -a (b15) 

—a 1 + f3]\j-i —a 

71 72 / 



where we define 
71 

73 

74 



72 



I 

2Ai' 
i i 
2Ai ^ 



giAx 
i i 
2At 



giAx 
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(B16) 



The main cause for artificial backrefiection in presence 
of the above boundary conditions comes from the approx- 
imate nature of the finite-difference evaluations (|B12[) 
and (jB13|) . Clearly, these approximations become bet- 
ter with decreasing grid spacing Ax, which means that 
a reduction of the grid spacing should lead to a more ef- 
ficient absorption of the outgoing flux. In practice, we 
find for grid spacings of the order of Ax — A/30 (with 
A = 27r/fc the wavelength of the condensate) that the 
relative amplitude of artificial backreflections from the 
grid boundaries is below 1%. We note that the amount 
of backrefiection that is accumulated during the numeri- 
cal propagation process would, at the same value of the 
grid spacing Ax, be considerably larger if the asymmet- 
ric version (jBlip of the one-way wave equation was used 
instead of Eq. (IBM]) . 
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